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Abstract 

A driven system of three species of particle diffusing on a ring is studied 
in detail. The dynamics is local and conserves the three densities. A simple 
argument suggesting that the model should phase separate and break the 
translational symmetry is given. We show that for the special case where the 
three densities are equal the model obeys detailed balance and the steady- 
state distribution is governed by a Hamiltonian with asymmetric long-range 
interactions. This provides an explicit demonstration of a simple mechanism 
for breaking of ergodicity in one dimension. The steady state of finite-size 
systems is studied using a generalized matrix product ansatz. The coarsening 
process leading to phase separation is studied numerically and in a mean-field 
model. The system exhibits slow dynamics due to trapping in metastable 
states whose number is exponentially large in the system size. The typical 
domain size is shown to grow logarithmically in time. Generalizations to a 
larger number of species are discussed. 

PACS numbers: 02.50.Ey; 05.20.-y; 64.75. +g 
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I. INTRODUCTION 



Collective phenomena in systems far from thermal equilibrium have been of considerable 
interest in recent years 0. Unlike systems in thermal equilibrium where the Gibbs picture 
provides a theoretical framework within which such phenomena can be studied, here no such 
framework exists and one has to resort to studies of specific models in order to gain some 
understanding of the phenomena involved. 

One class of such models is driven diffusive systems (DDS) |2|,|3|]. Driven by an external 
field these systems do not generically obey detailed balance so that the steady state has 
non- vanishing currents. Theoretical studies of DDS have revealed basic differences between 
systems in thermal equilibrium and systems far from thermal equilibrium. For example, it 
is well known that one dimensional (Id) systems in thermal equilibrium with short-range 
interactions do not exhibit phenomena such as phase transitions, spontaneous symmetry 
breaking (SSB) and phase separation (except in the limit of zero temperature or in the 
context of long-range interactions) 0]. In contrast, some examples of noisy Id DDS with 
local dynamics have been found to exhibit such phenomena. 

One example of a noisy system which exhibits SSB in Id is the asymmetric exclusion 
model of two types of charge studied in |5|||. In this model, two types of charge are biased 
to move in opposite directions on a Id lattice with open ends. The charges interact via a 
hard-core interaction, and are injected at one end of the lattice and ejected at the other end. 
This model is symmetric under the combined operations of charge conjugation and parity 
(PC symmetry). However, this symmetry is broken in the steady state, where the currents 
of the two charges are not equal. The reason for symmetry breaking in this model lies to 
some extent in the open boundaries. Other examples of models in which there is SSB in Id 
have also been found in the context of cellular automata @ and surface growth f|||. In 
the latter, SSB was due to the fact that one of the rates for a local dynamical move in the 
models is zero. Once this zero rate changes to a non-zero rate SSB disappears. 

A closely related problem to spontaneous symmetry breaking, is that of phase separation 
in Id noisy systems. This has been observed in driven diffusive models with inhomogeneities, 
such as defect sites [|TUJ or particles |TT| . In these models it has been found that macroscopic 



regions of high densities are formed near the defect, much like a high density of cars behind 
a slow car in a traffic jam [I2",I3|. Here the phase separation is triggered by the defects. 



It is of interest to study whether phase separation can occur in Id noisy homogeneous 
systems such as on a ring geometry with no defects, where all possible local transition 
rates which are consistent with the symmetry and conservation laws of the model are non 
vanishing. Recently, Lahiri and Ramaswamy have introduced a lattice model in the context 
of sedimenting colloidal crystals, where phase separation is found to take place without any 



inhomogeneities |L4]]. In this model, there are two rings coupled to each other and particles 
on each ring undergo an asymmetric exclusion process. The hopping rate between sites % 
and % + 1 on each ring depends on the occupation at the ith site on the other ring. However, 
this model is studied mainly using Monte Carlo simulations and no analytical results are 
available so far. 

In a recent Letter fl5| we introduced a simple three-species driven diffusive model ex- 
hibiting phase separation and spontaneous breaking of the translational symmetry on a 
ring. In the model nearest-neighbor particles exchange with given rates and the numbers 
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of each species are conserved under the dynamics. The rates of all local dynamical moves 
that obey the conservation laws are non zero. An argument indicating that generically the 
system phase separates, thus breaking the translational invariance, was given for the case 
when none of the species of particle has zero density. In the special case of equal number 
of particles of each type, it was shown that the local dynamics obeys detailed balance with 
respect to a long-range asymmetric (chiral) Hamiltonian. In this special case, using the 
Hamiltonian, we have found the steady state of the model exactly and have been able to 
prove the existence of phase separation analytically. 

The existence of a Hamiltonian for this special case is of interest in the light of spec- 
ulation that non-equilibrium systems exhibiting generic long-range correlations might be 
described by effective Hamiltonians containing long-range interactions [|16|Jl7| . Here we ex- 
plicitly demonstrate that for the special case where the three densities are equal the model 
is exactly described by a long-range asymmetric Hamiltonian. The model not only has long- 
range correlations but has generic long-range order. The mechanism found in this study 
suggests that systems with dynamical rules defined completely locally and a priori with- 
out respect to any Hamiltonian, may have a steady state where the configuration space is 
sampled according to a measure that is intrinsically global. The Hamiltonian also allows 
us to identify the analog of a temperature in the microscopic dynamics as related to the 
drive of the system; for zero drive, that is symmetric diffusion of the particles, the effective 
temperature is infinite and phase separation is lost. 

We note that a related but distinct three-species model has recently been introduced by 
Arndt et al. This model also exhibits phase separation 



In the present work we analyze in detail the M = 3 species model which was introduced 
in [15 1 and then generalize it to larger M. We provide the complete proof of phase separation 
which follows from the exact calculation of the partition sum in the thermodynamic limit. We 
also provide numerical evidence of phase separation in the general case where the densities 
of the three particles are not equal. 

In order to study the coarsening process Monte-Carlo simulations are performed. How- 
ever, simulation of the microscopic model is hampered by slow dynamics which makes it 
difficult to access the scaling regime. The system becomes trapped in metastable states 
comprising several domains of each type of particle. The number of metastable states is 
exponentially large in the system size. The lifetimes of the metastable states increase ex- 
ponentially with the average domain size as the fully phase separated state is approached. 
Thus the model provides an example of slow dynamics in a system without any quenched 
disorder fl9 |. 

To ameliorate the difficulty of numerically studying such slow dynamics we employ a toy 
model wherein it is the domains that are updated rather than the individual particles. This 
allows the long-time scaling behavior of the domain size to be investigated and to confirm 
a logarithmic growth of the average domain size with time. The toy model also affords a 
mean-field solution for the long-time dynamical behavior, that again confirms the scaling 
behavior. 

Returning to the case of equal numbers of particles of different species it is of interest to 
investigate the steady-state behavior in finite-size systems. We have found it convenient to 
do this by employing a matrix product technique previously used to solve the steady state 
of asymmetric exclusion processes ^D|. However, in the case of three species the simplest 
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form of this technique Jll ,|20, 12] is applicable only to a limited class of systems |21 1 . For the 



present model we generalize the matrix product to a product of rank 6 tensors and write 
the steady state by taking an appropriate contraction. The partition sum and steady-state 
correlation functions can be conveniently computed numerically using this tensor product 
ansatz. 



The paper is organized as follows: in section II we define the model introduced in [15 
and we present an argument which indicates that the system should phase separate as long 
as none of the species of particles has zero density. In section III we study the special 
case where the model satisfies detailed balance and explicitly write down the steady-state 
weight for the three-species model. The existence of phase separation in the model for any 
non-infinite temperature is proved analytically by calculating some bounds on the two-point 
correlation functions. Section IV contains numerical evidence for phase separation in the 
general case where the densities of the three species of particles are not equal. The toy 
model, which facilitates efficient Monte Carlo simulations, is used to study the dynamics of 
phase separation. A mean-field analysis of this toy model is presented, the details being left 
to Appendix A. In Section V we present results for finite systems obtained via the tensor 
product ansatz. Using these results we study finite-size scaling in the system. In section VI 
we address phase separation in systems with more than three species of particle and a proof 
of detailed balance for special cases is given in Appendix B. We conclude in section VII and 
discuss some open questions. 



II. DEFINITION OF THE MODEL 



We start by defining a three-species model which exhibits phase separation in Id. Con- 
sider a one-dimensional, ring-like (periodic) lattice of length N where each site is occupied 
by one of the three types of particles, A, B, or C . The model evolves under a random 
sequential update procedure which is defined as follows: at each time step two neighboring 
sites are chosen randomly and the particles at these sites are exchanged according to the 
following rates 

AB^=± BA 

BC^± CB (1) 
CA <0=± AC. 



The particles thus diffuse asymmetrically around the ring. The dynamics conserves the 
number of particles, Na, Nb and Nc of the three species. 

The q = 1 case is special. Here the diffusion is symmetric and every local exchange of 
particles takes place with the same rate as the reverse move. The system thus obeys detailed 
balance reaching a steady state in which all microscopic configurations (compatible with the 
number of particles Na, Nb and Nc) are equally probable. This state is homogeneous, 
and no phase separation takes place. We now present a simple argument suggesting that 
for q ^ 1 the steady state of the system is not homogeneous in the thermodynamic limit. 
For simplicity the case q < 1 is examined. As a result of the bias in the exchange rates 
an A particle prefers to move to the left inside a B domain and to the right inside a C 
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domain. Similarly the motion of B and C particles in foreign domains is biased. Consider 
the dynamics starting from a random initial configuration. The configuration is composed of 
a random sequence of domains of A, B, and C particles. Due to the bias a local configuration 
in which an A domain is placed to the right of a B domain is unstable and the two domains 
exchange places on a relatively short time scale which is linear in the domain size. Similarly, 
AC and CB domains are unstable too. On the other hand AB, BC and CA configurations 
are stable and long lived. Thus after a relatively short time the system reaches a state of 
the type . . . AAABBCC AABBBCCC ... in which A, B and C domains are located to the 
right of C, A and B domains, respectively. The evolution of this state takes place via a slow 
diffusion process in which, for example, the time scale for an A particle to cross an adjacent 
B domain is q~ l , where I is the size of the B domain. The system therefore coarsens and 
the average domain size increases with time as lnt/| lng] |[22| . Eventually the system phase 
separates into three domains of the three species of the form A. . . AB . . . BC . . . C. 

In a finite system the phase separated state may further evolve and become disordered 
due to fluctuations. However, the time scale for this to happen grows exponentially with 
the system size. For example it would take a time of order of q- mm { N B,N c } f or ^ e j± 
domain in the totally phase separated state to break up into smaller domains. Hence in 
the thermodynamic limit, this time scale diverges and the phase separated state remains 
stable provided the density of each species is non-zero. Note that there are always small 
fluctuations about a totally phase separated state. However, these fluctuations affect the 
densities only near the domain boundaries. They result in a finite width for the domain 
walls. The fact that any phase separated state is stable for a time exponentially long in the 
system size amounts to a breaking of the translational symmetry. 

Since the exchange rates are asymmetric, the system generically supports a particle 
current in the steady state. To see this, consider the A domain in the phase separated state. 
An A particle near the . . . AB . . . boundary can traverse the entire B domain to the right 
with an effective rate proportional to q Ne . Once it crosses the B domain it will move through 
the C domain with rate 1 — q. Similarly an A particle near the . . . CA . . . boundary can 
traverse the entire C domain to the left with a rate proportional to q Nc . Once the domain is 
crossed it moves through the B domain with rate 1 — q. Hence the net A particle current is 
of the order of q Ne — q Nc . Since this current is exponentially small in system size, it vanishes 
in the thermodynamic limit. For the case of Na = Nb = Nc, this argument suggests that 
the current is strictly zero for any N. In sections III and V we study this case in detail. 

The arguments presented above suggesting phase separation for q < 1 may be easily 
extended to q > 1. In this case, however, the phase separated state is BAC rather than ABC. 
This may be seen by noting that the dynamical rules are invariant under the transformation 
q — > 1/q together with A <-> B. 

III. SPECIAL CASE N A = N B = N c 

In this section we show that the dynamics (|1|), for the special case N A = Nb = Nc, sat- 
isfies detailed balance. The corresponding Hamiltonian, which determines the steady-state 
distribution, is found to have long-range asymmetric interactions. Using this Hamiltonian, 
we calculate analytically the partition sum and bounds on the correlation functions in the 
thermodynamic limit. These are then used to prove the existence of phase separation in the 
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model. Later, in section [V| we study finite systems for this case and the approach to the 
thermodynamic limit. 



A. Detailed Balance 

The general argument presented in the previous section suggests that for the special case 
Na = N B = N c , the steady state carries no current for any system size. We demonstrate 
this explicitly by showing that the local dynamics of the model satisfies detailed balance 
with respect to a long-range asymmetric Hamiltonian 7i. 

We define the occupation variables Ai, Bi and Ci as follows: 

^ j 1 if site % is occupied by an A particle 
1 1 otherwise. 

The variables Bi and Ci are defined similarly. Clearly the relation Ai + B ri + Cj = 1 is 
satisfied. A microscopic configuration is thus described by a set {Aj} = {Ai, Bi, Ci}. Using 
these variables, we will show that the Hamiltonian 7i and the steady-state distribution Wn 
corresponding to the dynamics ([T]) for the case Na = Nb = Nc = N/3 are given by 

N-l N 

n{{Xi}) = E E [CiB, - QAj + BiAj] , (3) 

i=l j=i+l 



W N ({Xi}) = Z N l q H ^ . (4) 

Here Zn is the partition sum given by J2Q n ^ Xi ^\ where the sum is over all configurations 
in which Na = Nb = Nc- Note that the Hamiltonian TC does not determine the dynamics 
of the system, it just governs the steady-state distribution as given in ([]) and (01). Eq. 
(|) suggests that q serves as a temperature variable with kT = — 1/lng. Thus, q — *■ 1 
is the infinite-temperature limit. The Hamiltonian (|3|) is written in a form which is not 
manifestly translationally invariant. However, careful examination reveals that when the 
relation Na = Nb = Nc is taken into account, the Hamiltonian as given by @ is indeed 
translationally invariant (see Appendix B). Therefore site 1 may be chosen arbitrarily. An 
expression for 7i which is manifestly translationally invariant will be derived at the end of 
this section. 
Note that 

JV-l N 

E E (C^i + A A) = (N/3) 2 (5) 

i=l j=i+l 

since the LHS yields the number of CA (and AC) pairs in the system. Using this relation 
the Hamiltonian may also be written in a form where the cyclic symmetry is more apparent: 

N-l N 

n({Xi}) = E E [CiBj + AiCj + - (N/3) 2 . (6) 

i=l j=i+l 
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The proof of Eqs. (|3],|]) is straightforward. This is done by considering a nearest-neighbor 
particle exchange and verifying that detailed balance is satisfied with respect to (§). We start 
by considering nearest-neighbor sites in the interior of the lattice, namely pairs other than 
(1,N). For example consider the exchange AB — > BA taking place at two adjacent sites k 
and k + 1, where k ^ N . This exchange results in the contribution of one more BiAj term in 
TC and hence the energy of the resulting configuration is higher by 1. It is easy to see using 
Eq. (f|) that qWjsi(- ■ ■ AB . . .) = Wn{- ■ ■ BA . . .), as required by detailed balance. Similar 
relations are easily derived for exchange of BC and CA pairs. Now consider an exchange 
taking place between sites N and 1, say CA — > AC. According to (|3p this exchange costs an 
energy of 2 Nb — Na — Nc + 1 • Therefore the exchange satisfies the detailed balance condition 
qW N (A . . . C) = W N (C . . . A) only when 2N B = N A + N C - Similarly by considering the 
exchanges AB — > BA and BC — > CB, one deduces that the detailed balance condition is 
satisfied for any exchange at sites A" and 1 as long as Na = Nb = Nc- In Appendix B 
we consider the most general nearest-neighbor exchange rates for M species and arbitrary 
densities and derive conditions (|B7|) for exchange rates which satisfy detailed balance. 

To write 7i in a manifestly translationally invariant form we define 7ii ({Xi}) as the 
Hamiltonian in which site io is the origin. Namely, 

N+i -2 N+io-1 

Wfc(W)= E E [CiBj - OA, + B^] , (7) 

i=i j=i+l 

where the summation over i and j is modulo N. Summing (0) over all io and dividing by 
N, one obtains, 

N N-l u 

K^Xi}) = E E (! - ^)(CiB i+k - C t A l+k + B^) (8) 

i=l k=l iV 
N N-l u 

= E E C 1 - M^ CiBi+k + AiCi+k + BiAi+k) ~ (iV/3)2 ' (9) 

i=l k=l ^ 

where in the summation the value of the site index (i + k) is modulo N. In the Hamiltonian 
(|) the interaction is linear in the distance between the particles, and thus is long-ranged. 
The distance is measured in a preferred direction from site i to site i + k. Moreover it is 
asymmetric in the sense that 7i is not invariant under the parity operation. It is also related 



to chiral Hamiltonians [23 



B. Ground States and Metastable States 

Before proceeding further to evaluate the partition sum and some correlation functions 
associated with the Hamiltonian (Q) let us make a few observations. The ground state of the 
Hamiltonian is given by the fully separated state A. . . AB . . . BC . . . C and its translationally 
related states. The degeneracy of the ground state is thus A^ and its energy is zero. A simple 
way of evaluating the energy of an arbitrary configuration is obtained by noting that nearest- 
neighbor (nn) exchanges AB — > BA, BC — > CB and CA — > AC cost one unit of energy each 
while the reverse exchanges result in an energy gain of one unit. The energy of an arbitrary 
configuration may thus be evaluated by starting with the ground state and performing nn 
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exchanges until the configuration is reached, keeping track of the energy changes at each step 
of the way. The highest energy is iV 2 /9 and it corresponds to the totally phase separated 
configuration A . . . AC . . . CB . . . B and its N translations. 

In considering the excited states of the Hamiltonian (|3]) we note that the model exhibits 
a set of metastable states which correspond to local minima of the energy: any exchange 
of nn particles results in an increase of the energy. In these states no BA, CB and AC 
nn pairs exist; only AB, BC and CA nn pairs may be found in addition to AA, BB, and 
CC. Any metastable state is thus composed of a sequence of domains in which A, B and C 
domains follow C, A and B domains, respectively. Therefore each metastable state has an 
equal number of domains, s, of each type with s — 1, ...,N/3. The s = 1 case corresponds 
to the ground state while s = N/3 corresponds to the ABC ABC . . . ABC state, composed 
of a total of iV domains each of length 1. (The total number of domains in an s-state is 3s.) 

For calculating the free energy and some correlation functions corresponding to the 
Hamiltonian (|3|) we find it useful to first derive some bounds for the number W(s) of s- 
states and their energies. In the following such bounds are presented. They are then used, 
in the next section, to evaluate the free energy and correlation functions of the model. 

To obtain a bound for jV(s) we note that the number of ways of dividing N/3 A particles 
into s domains is [ N ^j • The number of ways of combining s divisions of each of the three 

types of particles is clearly ( JV /^ 1 ) ' There are at most iV ways of placing this string of 
domains on a lattice to obtain a metastable state (the number of ways need not be equal to 
N since the string may possess some translational symmetry). One therefore has 



'N/3 - r 
. s-1 . 



< Af(s) < N 



'N/3 - 1 N 
. s-1 . 



1 3 



(10) 



Thus, the total number of metastable states is exponential in N. 

We now consider the energy of the metastable states. It is easy to convince oneself that 
among all s-states, none has energy lower than the following configuration, 



A...AB ...BC ... C ABC ABC . . . ABC 



(11) 



where the 3(s — 1) rightmost domains are of size 1 and the three leftmost domains are of 
size {N/3 — s + 1) each. The energy of this state, E s satisfies the following recursion relation 



E s = £ s _x + N/3 - s 



(12) 



with Ei = 0. To see this one notes that the s-state may be created from the (s — l)-state 
by first moving a B particle from the leftmost B domain across (N/3 — s) C particles to 
the right and then moving an A particle from the leftmost A domain to the right across the 
adjacent B and C domains. The energy cost of these moves is (N/3 — s), yielding (0). The 
recursion relation (|1^), together with E\ = 0, is then readily solved to give 



1) — 

1 3 



s(s-l) 



(13) 



The energy of all metastable s-states is larger or equal to E s as given by Eq. (|T3"D . In the 
following Section we use the bounds (|T0| ) and ([13]) to calculate the partition sum and some 
correlation functions corresponding to the Hamiltonian (W). 



S 



C. Partition Sum 



In this section we prove that, in the large N limit and for all q < 1, the partition sum is 
given by 

Z N = iV/[(g)oo] 3 , (14) 

where 

{q) 00 =Vmfl-q){l -<?)...(!-<?). (15) 



The partition sum for q > 1 may be obtained by replacing q by \jq in fll~4]). Note that 
the partition sum is linear and not exponential in N, meaning that the free energy is not 
extensive. This is a result of the long-range interaction in the Hamiltonian and the fact that 
the energy excitations are localized near the domain boundaries, as will be shown in the 
following. 

For q close to 1, (q)oo has an essential singularity, 

(g)oo = e-i^ 2 / 6 ^ 1 ^] . (16) 

This suggests that extensivity of the free energy could be restored in the double limit q — ► 1 
and N — > oo with iVlng finite. This scaling behavior is only suggestive since expression 
(|14]) may not be valid in this limit. For q = 1 all configurations with = Nb = Nc are 
equally probable, so that the partition sum is given by Zn = (jvy 3 ) (^n/s) which g° es like 
3^ for large N. 

It is instructive to first present the proof of flI3] ) for small values of q. This proof will 
then serve as the basis for the proof for any q < 1. 

We start by noting that in calculating the partition sum fll4|) , configurations with energy 
larger then aN (a > 0) may be neglected in the thermodynamic limit. For simplicity we first 
demonstrate this for q < (l/S) 1 ^, although later we show it for any q < 1. The contribution 
to the partition sum from these energy states, Z m>aN , is given by 

N 2 /9 

Z m >aN= E D(rn)q m , (17) 

m=aN+l 

where D(m) is the number of configurations of energy m. Clearly, the number of possible 
configurations in the system is bounded crudely from above by 3^. This bound implies 

N 2 /9 

Z m> aN< E SV"- (18) 
m=aN+l 

Thus, for q < (l/S) 1 ^, the contribution to the partition sum arising from energies larger 
than aN is exponentially small in iV and may be neglected in the thermodynamic limit. 

The calculation of Z^ is thus reduced to calculating a truncated partition sum in which 
only energies up to aN are summed over. To proceed we consider a < 1/3 and take 
into account configurations with energy less than iV/3 — 1. This simplifies the calculations 
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FIG. 1. The I = 5 ground state for an N = 21 system. 



considerably since all configurations with energy m < N/3 — 1 may be decomposed into 
N disjoint sets of states, each corresponding to a unique ground state (see previous section 
for a discussion of ground states). We label the sets by I = 1...N, the position of the 
rightmost A particle in the A domain of the ground state belonging to the set (see Fig. [I] 
for an example of an I = 5 ground state). Each state in a specific set can be obtained from 
the corresponding ground state by exchanging nearest neighbors so that the energy always 
increases along the intermediate states. Note that this is correct only if excitations of energy 
less than N/3 — 1 are considered. This is because not all higher energy states can be reached 
by uphill steps from a ground state. 

Thus, using translational invariance, the partition sum can be written as 

Z N = NZ N + e-°W , (19) 

where Zn is the truncated partition sum of one of the N sets of configurations. 

We proceed to calculate Zn- This is done by considering all the possible energy exci- 
tations with energy less than N/3 — 1 above one ground state. Consider a specific domain 
boundary, say AB. Excitations of energy m at this boundary can be created by moving 
one or more A particles into the B domain (this is equivalent to moving B particles into 
the A domain). An A particle moving into the B domain is considered as a walker. The 
excitation energy increases linearly with the distance the walker has moved. Thus, in this 
picture an excitation of energy m is created by 1 < j < m walkers, traveling a total distance 
m. The number of excitations of energy m is then given by the number of ways, P(m), 
of partitioning an integer m into the sum of a sequence of non-increasing positive integers. 
Taking into account excitations at all three boundaries, an excitation of energy m in the 
system is created by three independent excitations of energy mi, 777,2 and 772,3 at the different 
domain boundaries such that mi + 7712 + 7713 = m. The number of excitations of this form 
is just given by P(mi)P(m 2 )P{m 3 ). Thus, Z N is given by 

JV/3-2 m 

Z N = ^2 <? m E P{mi)P(m 2 )P (20) 

771=0 771i=0 

Taking the thermodynamic limit we obtain 

00 

lim Z Ar = (^g m P(m)) 3 . (21) 

777 =0 

Using a well known result from number theory, attributed to Euler, for the generating 



function of Pirn) [24 
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E ? m ^M = ( 22 ) 



m=0 



and using and (|19D, Eq.(|14]) is obtained. 

So far we have proved that for g < (1/3) 3 , Eq.(14) is exact in the thermodynamic limit. 
We now extend these results for any q < 1. First we have to show that the states ignored in 
the previous calculation for q < (1/3) 3 may be ignored for all q < 1. To do this we calculate 
upper and lower bounds on Zn and show they converge for large enough N. 

For this we have to consider the entire energy spectrum of the Hamiltonian. Any config- 
uration of the system which is neither a ground state nor a metastable state can be obtained 
from at least one ground state (s = 1) or a metastable state (s > 1) as follows: starting 
from this s-state exchange nearest neighbors such that the energy always increases along 
the path until the configuration is reached. In what follows it is demonstrated that none of 
the configurations which can be obtained from s-states, with s > 1, by the above procedure 
of particle exchange, contributes to the partition sum in the thermodynamic limit. 

An upper bound on the partition sum may be calculated as follows: using the same 
steps of derivation used for computing Zn, it is straightforward to show that the con- 
tribution to the partition sum from an s-state and associated configurations is at most 
^(s-i)iV/3-s(s-i)/2|-^^j-3s_ rj^g p re f ac tor g( s_1 ) 7V / 3_s ( s_1 )/ 2 arises from the minimum energy 

(|13|) of this metastable state. Therefore by considering the contributions from all the s-states 
and using ( |TUD the following bound is found 

Z N < N/[(qU 3 + 1 )? (s - 1)Ar/3 - s(s - 1)/2 [(?oo)]- 3s ■ (23) 

The second term on the RHS represents the contribution from excitations around the 
metastable states. Replacing g( s_1 ) iV / 3 ~ s ( s_1 )/ 2 by an upper bound g^" 1 )^/ 6 one can sum 
the binomial series. The resulting expression is exponentially small in N for any q < 1. 

A lower bound on can be calculated by neglecting configurations with energy greater 
than iV/3 — 1 as follows, 

N/3-2 m 

Z N >N <T E P(m 1 )P(m 2 )P(m 3 )5 mi 

+m2+m,3,m 

(24) 

m=0 m,i=0 

oo m 

= N/Kq)^] 3 -N V m E P{m 1 )P{m 2 )P{m 3 )5 mi+m2+m3 , m (25) 

m=JV/3~l mi=0 
oo 

> N/Kq)^ 3 - N £ q m (mP(m) f . (26) 

m=N/3-l 



The asymptotic behavior of P{m) |24| is given by 

P( m ) ~ -J— exp (vr(2/3) 1/2 m 1/2 ) . (27) 
4m v 3 

Thus, for large N the lower bound (|26| ) converges to ([14]) as does the upper bound (|23|). 
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D. Correlation Functions 

Whether or not a system has long-range order in the steady state can be found by 
studying the decay of two-point density correlation functions. For example the probability 
of finding an A particle at site i and a B particle at site j is, 

(ABj) = ^ £ AiB s q n « x »» , (28) 

where the summation is over all configurations {Xk} in which Na = Nb = Nc- Due 
to symmetry many of the correlation functions will be the same, for example (AiAj) = 
(BiBj) = (CiCj). A sufficient condition for the existence of phase separation is 

lim lim ((AM - (A 1 )(A r )) > 0. (29) 

Since (Ai) = 1/3 we wish to show that lim^oo limjv->.oo(^4iA-) > 1/9. In fact we will show 
below that for any given r and for sufficiently large N, 

(A 1 A r ) = 1/3 - 0(r/N) . (30) 

This result not only demonstrates that there is phase separation, but also that each of the 
domains is pure. Namely the probability of finding a particle a large distance inside a domain 
of particles of another type is vanishingly small in the thermodynamic limit. 

To prove Eq. (|3"CD, we use the relation (AiA r ) = 1/3 — (A\B r ) — (AiC r ), and show 
that the correlation function (AiB r ) is of 0(r/N) and (AiC r ) is of Q(l/N). Here we show 
only the proof for (A\B r ), since the proof of (A\C r ) is similar. We also restrict ourselves to 
r < N/3, which is sufficient for proving Eq. (|30|). 

We have already seen that the contribution to the partition sum from the metastable 
states and excitations above them are exponentially small in the system size and hence may 
be neglected. Therefore, for calculating the correlation function it is sufficient to consider 
the N ground states and excitations above them which may be reached by moves which 
only increase the energy. As we have seen, these states form iV disjoint sets of states, each 
associated with one of the ground states. Using this we now show that (AiB r ) = 0(r/N). 
For this purpose we use a restricted partition sum Z s , which is defined as the partition sum 
Zn calculated with the constraint that one of the walkers, say of type A, has traveled at 
least distance s. It is given as iV — > oo by 

oo m 

Zs=Y. <? m E P s (m 1 )P(m 2 )P(m 3 )5 m (31) 

m=0 m,i=0 

Here P s (m) is the number of partitions of integer m with the constraint that in all the 
partitions the integer s occurs at least once. Noting that P s (m) = P(m — s) it is easy to 
show that 

Z s = q s Z , (32) 

where Z = liniAr^oo Z N . 
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1 r I N 

FIG. 2. A ground state where A\ = 1, B r = 0. In order to have B r = 1, it is necessary for a B 
particle to travel a distance of I — r + 1 into the A domain. 

We now proceed to derive a bound for (A%B r ). Recall that I is the position of the 
rightmost A particle in the A domain in the ground state labeled I. If we define [A\B r )i as 
the correlation function calculated within the set of states labeled I, we can write 

1 N 

= -^(AiBji , (33) 
iy i=i 

up to exponentially small corrections in the system size. For convenience we break the 
summation over I into 4 sums according to the values of A\ and B r in the ground state. 
These 4 parts correspond to (I) ground states where A\ = l,B r = 1; (II) ground states 
where Ai = l,B r = 0; (III) ground states where A 1 = 0,B r = 1 and (IV) ground states 
where A\ = 0, B r = 0. We now consider each of them in detail and give an upper bound for 
(A\B r )i in each case. 

(I) Ground states where A± = 1, B r = 1: In this case the site 1 is inside the A domain and 
site r is inside the B domain. Since we consider only r < N/3, these states correspond to 
the / ground states with 1 < I < r. Using the fact that (A\B r )i < 1 one finds 

Y d {A 1 B r ) l < r - 1 . (34) 
i=i 

(II) Ground states where A\ = 1, B r = 0: in principle site r might be either inside the A 
domain or inside the C domain. However, since site 1 is in the A domain and we consider 
only r < N/3, site r must be in the A domain. The ground states I for which this takes place 
satisfy r < I < N/3. Clearly, only the excited states where B r — 1 contribute to {A\B r )i. 
In such excited states one of the B walkers travels at least a distance I — r + 1 into the A 
domain (see Fig. §). F° r this case we can give the upper bound (A\B r )\ < J2t=i- r +i %sl 
From Eq. ©, Z,/Z = q s . Hence, 

N/3 N/3 oo 

E<^i2*r>i<E E <7 S - (35) 

l=r l=r s=l— r+1 

(III) Ground states where A\ = 0, B r = 1: again, since r < iV/3 site 1 has to be inside the 
B domain. The values of / satisfying this condition are in the range 2N/3 + r < I < N. 
In this case only excited states in which one of the A walkers travels at least a distance 
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N — I + 1 into the B domain will contribute to {A\B r )i. Hence we can use the upper bound 
{AiB r )i < Y^T=n-i+i^s/ in this case. Therefore, 

N N oo 

£ (A\B r )i < Y, E I s - (36) 

l=2N/3+r l=2N/3+r s=N-l+l 

(IV) Ground states where A\ — 0, B r — 0: there are three possibilities here, (a) site 1 is 
inside the C domain and site r is inside the A domain (iV/3 < I < N/3 + r), (b) both the 
sites 1 and r are inside the C domain (r + iV/3 < I < 2N/3), (c) site 1 inside the B domain 
and site r inside the C domain (2N/3 < I < 2N/3 + r). Since all these are consistent with 
r < N/3, all these cases can occur. It can be shown that the minimal energy needed to 
create an excited state where A\ = 1 and B r = 1 is e a = 21 — r — N/3 — 1 for the case (a), 
6f) = N/3 + r — 3 for the case (b) and e c = 5iV/3 — 21 + r — 1 for the case (c). The resulting 
expression for the bound is 

2AT/3+r-l N/3+r-l oo 2AT/3 oo 27V/3+r-l oo 

E (^A.),< E E? s + E E? s + E E? s - (37) 

/=AT/3+l Z=/V/3+l s= e f £=/V/3+r s=e& /=2Af/3+l s=e c 



The summations on the RHS of Eqs. (pq - pTD can be carried out explicitly. To leading order, 
the summations gives q/(l — q) 2 for each of Eqs. ( |3"5| ) and The summation on the RHS 
of Eq. (^) vanishes exponentially in the thermodynamic limit. Using Eqs. (|33| - |36D , we get 
the following expression for the upper bound on (AiB r ) 



2? , _-0(iV) 



1+ (^ + e 



(38) 



Therefore (AiS r ) = 0(r/N). Similarly one can show that (AiC r ) = 0(1/ N). Thus for all 
g < 1, (AiAf) = 1/3 — 0(r/N), proving the existence of a complete phase separation. 



IV. COARSENING 
A. Monte Carlo Simulations 

We have demonstrated that in the thermodynamic limit the system is phase separated 
when Na = Nb = No- The general arguments given in section II indicate that when the 
global densities of the three species are non- vanishing and q ^ 1, the system phase separates, 
even when the three densities are not equal. The argument suggests that the typical time, tf, 
in which the system leaves a specific phase separated configuration increases exponentially 
with the system size. Thus, a phase separated state is stable in the thermodynamic limit. 
In the following we use Monte Carlo simulations to support these arguments. 

The time, tf, can be measured using the auto-correlation function defined as, 

1 N 

C W = 77 E((MQ)Mt)) + (B t {0)B t (t)) + (Ci(0)tf,(t)» , (39) 
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where Ai(t),Bi(t) and Ci(t) are the values of the occupation variables A^Bi and Cj at 
time t, and (...) denotes an average over histories of evolution. Clearly, c(0) = 1 while 
c(oo) = (N A /N) 2 + (N B /N) 2 + (N c /N) 2 , the value of the autocorrelation between two 
independent configurations. Thus, tf may be defined as the decay time of c(t) to c(oo) when 
at t — the system is totally phase separated. 

We have measured the time scale tf using Monte Carlo simulations for different system 
sizes for N A = iV^ = Nc and for Na 7^ Nb 7^ Nc for several q values. An example of such 
measurements for Na/N = 0.4, Nb/N = 0.35 and Nq/N = 0.25 is presented in Fig. |3|. In 
the plot tf is plotted versus system size for several values of q. It agrees with the exponential 
growth of tf with the system size suggested by the simple argument of section II. The same 
behavior seems to occur for all q ^ 1 and different choices of Na/N, Nb/N and Nc/N. 
Therefore we conclude that the Monte Carlo simulations support the claim that for any 
q 7^ 1 the system will phase separate into three domains in the thermodynamic limit, even 
when the number of particles of each species is not equal. In the thermodynamic limit the 
translational symmetry is spontaneously broken in this state. Due to the slow dynamics, 
which reflects escape from metastable states, Monte Carlo simulations could be performed 
only for relatively small system size (N ~ 100). In order to study the coarsening process for 
larger systems we employ, in the following, a toy model which mimics the dynamics of the 
model ([]]). The toy model can be conveniently simulated for systems larger by about two 
orders of magnitude. 

B. Toy Model 

We now construct a simple toy model which captures the essential physics of the coarsen- 
ing process in the model at large times and enables us to simulate systems much larger than 
those accessible by Monte Carlo simulation. Using the toy model we examine another char- 
acteristic scale of the system. Namely, the average domain size (/) as a function of the time t. 
The results support the simple argument leading to a domain growth law (I) ~ log t / 1 log g | . 
A mean-field version of the toy model is then solved analytically. 

We consider a system at time t such that the average domain size, (I), is much larger 
than the domain wall width. At these time scales, the domain walls can be taken as sharp 
and we may consider only events which modify the size of domains. This means that the 
dynamics of the system can be approximated by considering only the movement of particles 
between neighboring domains of the same species. Using this we represent a configuration by 
a sequence of domains of the form A 1 B 1 CiA 2 B 2 C 2 . • . A k ~B k Ck, where the ith domain of, 
say A, particles is represented by Aj, as shown in Fig. The exchange of particles between 
domains, say A« and Aj+i, takes place at a rate dictated by the size of the domains Bj and 
Cj which separate them. Since intermediate configurations of the form . . . Aj_iBjAjCj . . . 
rearrange on short time scales compared with the evolution between metastable states, only 
metastable configurations are considered in the toy model. Events in which a domain splits 
into two are ignored. 

Using these ideas we define the dynamics of the toy model as follows: at each time step 
two neighboring domains of the same species of particle are chosen randomly, say Aj and 
A i+1 . Let dj, bi and q denote the lengths of the domains Aj,Bj and Cj respectively. The 
length of the domain chosen is then modified by carrying out one of the following processes: 
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FIG. 3. The decay time tj measured for different system sizes for several q values. Here 
N A /N = 0.4, N B /N = 0.35 and N c /N = 0.25. The data is averaged over a 100 runs. 



1) at 



Oj+l —> CLi+l + 1 



2) a, 



fli+l a i+l — 1 



with rate q b 



with rate q c 



(40) 



where, as before, q < 1 is considered. 

If a; becomes zero, then delete the domain Aj from the list of domains, and merge Bj 
and Cj with Bj_i and Cj_i, respectively. Then for j > i, shift the indices of the domains 
from j to j — 1, so that K becomes K — 1. The rules for updating B and C domains can 
be obtained from (|4Tj| ) using cyclic permutations and a slight change of indices. 
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FIG. 4. A configuration of the toy model represented by a sequence of domains. 
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Note that the toy model is only relevant to the description of the coarsening dynamics. 
This is because here, once the system is left with three domains, it remains in that state. 

To simulate the toy model efficiently, an algorithm suitable for rare event dynamics must 
be used due to the small rate of events f2"5fl . We use an algorithm which is performed by 
repeating the following steps: 

1. List all possible events {n} and assign to them rates {r n } according to the rules of the 
model. 

2. Choose an event m with probability r m /R where R = J2 n r n- 

3. Advance time by t — > t + r, where r = l/r m . 

The algorithm would be equivalent to a usual Monte Carlo simulation, where 1 time step 
is equivalent to one Monte Carlo sweep, if in step 3, r would be drawn from a Poisson distri- 
bution R exp[—Rr]. However, here we make an approximation by using the deterministic 
choice t = l/r m . 

We have simulated the dynamics for lattices of size up to 9000. For simplicity we consider 
the case Na = = Nq- An example of a typical behavior of the average domain size as a 
function of t is shown in Fig. [5]. One can see that after an initial transient growth time the 
data fits very well with a log(i) behavior. Simulations for different q values indicate that, 

(0 = alogt/|logg| (41) 

with a ~ 2.6. The toy model enables one to verify the scaling behavior (|4"ID and estimate 
the constant a. This would be very difficult to do by simulation of the full model (U). 



C. Mean Field Solution of the Toy Model 

Here we present the solution of a mean-field version of the toy model based on ideas 
presented by Rutenberg and Bray p6 and Derrida et al. |27| in the study of the ordering 



dynamics in a one- dimensional scalar model. To construct the mean-field model we notice 
that since all steps in the toy model which involve exchange of particles between domains 
occur at a rate exponentially small in the size of domains, one can consider a model where 
dynamics occurs only in the vicinity of the smallest domains. The mean field approximation 
assumes that different domains are uncorrelated and does not distinguish between domains 
of different species. The second assumption relaxes the conservation of particles of each 
species. Thus, in contrast to the systems studied by and |27J we do not expect the mean 
field to become exact in the scaling limit. We define the mean- field model as follows: 

1. Pick one of the smallest domains D m ™. 

2. Pick 2 domains D\ and D2 randomly and treat them as the neighbors of -D m i n . 

3. Pick 3 more domains randomly say, D 3 , and D 5 . 

4. Eliminate -D m ; n , D2 and D3 from the system and change the length of the domains 
D 3 , D 4 and D 5 by 
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FIG. 5. Monte Carlo simulation results for the toy model for the average domain size, (I), vs. 
time ,t, for N = 9000 and q = 0.8. The data is averaged over 1760 runs. 



Id 3 Id 3 + Id 2 

ln r — > Inc. + In ■ 



(42) 



Steps 1 to 4 are performed simultaneously for all of the smallest domains in the system. 
Here lo t is the length of domain i, and Z m i n is the length of the smallest domain. Steps 1 to 
3 choose the smallest domain D min and its nearest neighbors (see Fig. g). Step 4 uses the 
fact that the dynamics occur in the model only in the vicinity of the smallest domain and 
eliminates the three domains D min , D\ and D 2 joining them appropriately with the other 
domains D 3 , _D 4 and D 5 . Note that this mean- field dynamics does not take into account the 
time taken for these events to happen. This will be done later when we derive the growth 
law of domains. 

To solve the mean-field model we follow the method used in p6|, 27 1 . Let ni(t) be the 
number of domains of size /, irrespective of the type of particle it consists of. Let / m i n (£) 
be the length of the smallest domain and M(t) = 3K(t) be the total number of domains at 
time t. We will assume that ni(t) has the following scaling form in the large-time limit. 
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FIG. 6. A configuration in the mean- field model after steps 1 to 3 have been performed for one 
of the smallest domains. 

«, = ^f(j^) ■ (43) 

^min ^min 

A solution of the model given in Appendix A yields (I) = (x)l m i n where (x) is given by 

3e 7/3 

^ = / °° dx x-V3 e -x e -i{x)/3 ■ ( 44 ) 



The growth law for Z mm (t) can be derived following p6[: after the elimination of the 
smallest domain, Z m ; n increases by 1. This happens at a rate q linin /(l), namely, the inverse 
time required by a typical domain to cross a distance i m i n (thus, causing the annihilation of 
Anin)- Using (I) = (x)l m i n we write 

f)l . drain 

L/l, min y (45) 

From this one can obtain the scaling form of the average domain size, 



Hq)\ 



t . (46) 



Note that in this Eq. (x) does not depend, according to the mean-field solution, on either 
t or q. One can see from fl46|) that for large /, (I) ~ (x) lnt/| lng| which was confirmed by 
the simulations of the toy model where a = (x) (see (f4l~D). A numerical evaluation of 
yields (x) ~ 3.72 as compared with a = 2.6 obtained from the toy model simulations. 



V. EXACT RESULTS FOR FINITE SYSTEMS 

In section III the partition function, Z^, and correlation function, (A\A r ) for finite r, 
have been calculated in the thermodynamic limit. It is also of interest to obtain results for 
finite systems for the study of finite-size effects and the approach to the thermodynamic 
limit. Recently a matrix ansatz method has been introduced to study one-dimensional 
non-equilibrium systems pOj. It has been shown that in certain three-species models the 



steady-state weight and correlation functions can be represented as a product of matrices 
|TT|,|5| JT2]j2T[| . In the ansatz a specific matrix is associated with each type of particle. Then 
the unnormalized probability of a certain configuration is obtained from a matrix product. 
The matrices corresponding to the different species of particles satisfy an algebra derived 
from the dynamics of the model. A scalar, i.e. the weight of a configuration or some 
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correlation function, is usually obtained by performing a trace over the product of matrices, 
or by multiplying both sides of the product of matrices by vectors. Generalizing this method 
to replace matrices by tensors |[28|| , we have been able to obtain recursion relations for the 
partition function and correlation function for finite systems for the special case Na = Nb = 
Nc- The recursion relations are then used to obtain the partition function and correlation 
function (AiA r ) for any r in small systems. The results are used to study the scaling of the 
correlation function near the critical point q — 1 (infinite temperature), where the typical 
domain wall width diverges. 



A. The Tensor Product Ansatz 

It is convenient to consider the unnormalized weights, f^({Xi}), defined through 

W N ({Xi}) = Z^f N {{Xi}) , (47) 

where Ww({Xi}) is the probability of being in configuration {Xi}. The partition sum Z N 
is given by 

Zn=Y. MTO) , (48) 

{Xi} 

where the sum is over all configurations with Na = Nb = Nc- 

We generalize the matrix ansatz and construct the steady-state weight, /jv({-Xj}), from a 
product of tensors each corresponding to a particle located in a specific place on the lattice. 
The contraction of the tensors yields a tensor which is then contracted with 'left' and 'right' 
tensors to generate a scalar. The three tensors which represent the different type of particle 
are defined as rank 6 tensors through the following tensor products of square matrices 

A = E<8>D®1 

B = 1 <g> E <g> D (49) 
C = D® 1® E . 

Here 1 is a unit matrix. The matrices D and E will be chosen in what follows to satisfy a 
commutation relation which will be dictated by the detailed balance condition. 

To define /jv({^Q}) we introduce the following notation: the contraction of two rank 6 
tensors O = Oi <E> 2 <8> 3 and "P = Pi <g> P 2 <8> P3, where Oi and Pi are square matrices, 
according to the rule O1P1 ® O2P2 <8> O3P3 is denoted by OV . The contraction of a rank 
6 tensor O with a 'left' rank 3 tensor (/C| = (K\\ ® (K 2 \ ® {K 3 \, where (Ki\ are transposed 
vectors, and a 'right' rank 3 tensor, \M) = \Mi) ® \M 2 ) ® I-M3), where |Mj) are vectors, 
defined through ( J Rr 1 |0 1 |Mi)(K 2 |0 2 |M 2 )(ir3|03|M3) is denoted by (/C| O \M). 

Using these definitions we write the steady-state weight of the system as 

N 

MiX,}) = (U\ UlAA + B.B + dC ]|V) , (50) 
1=1 

where A^, Bi and Cj are the occupation variables defined in (0). The expression states that 
a tensor A is present at place i in the tensor product if site % is occupied by an A particle, a 
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tensor B is present if site i is occupied by a B particle, and a tensor C is present if site i is 
occupied by a C particle. The action of the tensor product on (U\ and |V) (to be specified 
later) produces the scalar f^({Xi}). 

It is straightforward to show using detailed balance that a necessary condition for Eq. 
(p0|) to be the steady-state weight is that the following commutation relations are satisfied 
between A, B, and C : 

q AB = BA 

qBC= CB (51) 
q CA= AC . 

Using fl^Df ) one can verify that these commutation relations are satisfied provided that gDE = 
ED. This deformed commutator is of relevance in other stochastic systems p9| , p0| . A 
representation of the matrices D and E which satisfies this commutation relation can be 
obtained as follows: let {{n\} denote a basis set (n = 0, 1, . . . , N/3) forming a vector space. 
In this basis we choose the matrices so that, 

(n|E = (n\q n for any n 
(n|D — (n — 1\ for n > 1, 

while for n = 0, (0|D = 0. An explicit form for E and D is given by the following (N/3 + 
1) x (N/3 + 1) square matrices 

N/3 N/3 

E=^2l n \n)(n\ ; D = £ \n)(n - 1| . (53) 

n=0 n=l 

To obtain f^({Xi}), (U\ and |V) have to be specified. This should obviously be done so 
that fpf({Xi}) is non-zero if the ansatz is to give a non-trivial result. We consider a general 
tensor product which corresponds to some configuration. The product has N/3 tensors of 
each type A, B and C , which results in a tensor product of three matrix products. Using 
( |49|) it can be seen that each matrix product contains N/3 matrices of each type D, E and 
1. Since E and 1 are diagonal, while D acts to its left as a lowering matrix (see (|5"2"D), 
choosing (U\ = (N/3\ <8> (N/3\ ® (N/3\ = (N/3, N/3, N/3\, and |V) = |0,0,0) will give a 
non-zero fw({Xi}). This makes clear that the minimal size choice for the vector-spaces is 
N/3+1. Under this choice it is easy to see that in the ground states one has = q N ' 2 / 9 which 
corresponds to a ground state energy N 2 /9. However, the choice of (U\ and (V| is determined 
only up to some multiplicative factors. These factors may be used to shift the ground state 
energy of the system. For example choosing (U\ as before with |V) = q~ N ' 9 |0,0, 0) will 
shift the ground state energy to 0. In the following the factors are taken to be 1. 

Finally we would like to remark that usually when using the matrix ansatz for systems 
with periodic boundary conditions it is often convenient to use a trace of the matrix product 
PUpII . In this case this is not possible since the trace of our tensor product is always zero. 

B. Partition Sum 

In terms of the ansatz (|50|) the partition function Zj^ is given by 
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(52) 



,N N N, 



Zn = (j , - , - \(A + B + C)»\0, 0, 0) . (54) 

Note that we are using the canonical ensemble since all tensor products with unequal number 
of particles do not contribute to Zn- This is easily seen since in these cases there are always 
more than iV/3 D matrices acting on one of the vectors (N/3\. 

To obtain the partition function we derive a recursion relation for 



G l ^ k = (i,j,k\(A+B + C) l \0,0,0) . (55) 
One can see that Gn/3,n/z,n/z = Rewriting G\j k as 

G\ hk = (i.j, k\ A(A +B + Ct'lO, 0, 0) + k\B(A+B + C )'- 1 |o, 0, 0) 

+(i,j,k\ C(A+B + C ) | - 1 |0, 0, 0> , (56) 

and using relations Q4"5| ) and ( p2] ) the following recursion relation can be derived 

^ij,k = <fG l i,j-i,k + tfGijjk-i + Q k G l i-i,j,k ■ (57) 

The boundary conditions for this recursion relation is given by the no particle partition 
function G° • k = lifi = j = k = and is zero otherwise. 

For small systems (up to N = 21) for which the recursion relation is tractable analytically 
on Mathematica, we obtained the partition function Zn = Gj^, 3 N , 3 N , 3 as a polynomial in 
q. As expected the first N/3 — 2 terms of the polynomial match the first N/3 — 2 terms of 
the expansion of fll4]) up to a factor of q N / 9 due to the energy shift in the ground state. For 
larger A" we solve the recursion relation numerically. 

We note that (0) could have been derived directly from the definition of the partition 
function without recourse to the tensor ansatz. However, we believe the utility of the ansatz 
lies in the ease with which correlation functions can be manipulated and relations such as 
that of the next subsection derived. 

C. Correlation Functions 

The correlation function (A\A r ) is given in terms of the ansatz by 

{AiA) = (ff, f, f I A{ A + B + C) r - 2 A(A+B + C) N - r \0, 0, 0) _ 

Zn 



Using relations (|49|),(|5T|) and (|52| ) we obtain 

{AiA) = q ^ (f ; f ~ 2, f \{A + qB + C (ir%A +B + C )"-|0, 0, 0) (gg) 

Zn 

= q~ JH 8 -^, (60) 

ZjN 

where we define the object U(r)*j k through, 
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FIG. 7. The correlation function C = (AiAjy/z), obtained from the tensor ansatz, as a function 
of the scaled variable iVlng for N = 30,36,42, . . . ,84. The inset shows the same data plotted 
against q. 

U(r)l jik =(i,j,k\(A + qB + C /q) s - 2 (A+B + C ) N - S \0,0,0) . (61) 

A recursion relation for U(r)fj k can be obtained, similarly to the recursion relation for the 
partition function (|5"TD . Using (|4*9|) and ( (52]) gives 

H'-);,., = q l U(r)t-\ k + q j+1 U(r)^U + ^(r)^ . (62) 

The boundary conditions for the recursion relation are obtained by noting that U (r)? • fc = 

{i,j,k\ (A- + B + C) N ~ r \0, 0, 0), i.e., U(r)fj k = G^~ k . Using the same methods one can 
obtain recursion relations for all other correlation functions. 

The recursion relations are solved numerically for finite systems. This is done by first 
solving numerically for Gf~ k , and then using the result as boundary conditions for the 
recursion relation fl62]). Owing to (|60|) we are ultimately interested in U(r) r N , 3 N / 3 _ 2 n/3- 

Using these recursion relations we have calculated the correlation function (AiAn/2), 
which is a measure of the phase separation in the system, for various system sizes. When 
(AiA N / 2 ) is close to zero the system is phase separated. In the disordered case, q — 1, the 
value of the correlation function is (N — 3)/9(iV — 1), approaching 1/9 in the thermodynamic 
limit. The results are shown in the inset of Fig. [7] one can see that the system is phase 
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separated for small values of q while for q close to 1 the system is disordered. The range of 
q values for which the system is phase separated increases as the system size increases. The 
natural scaling variable near the critical point q = 1 is Nlnq, the ratio between the domain 
size A/3 and the domain wall width Jlq l dl/ J q l dl = l/|lng|. In Fig.(^) the correlation 
function is plotted as a function of the scaling variable. One can see that the data collapse 
improves as the system size increases. This scaling variable was also suggested by the form 
of the partition function (|TJ) in section III C. 



In the following we discuss possible generalization of the model to M > 3 species. To 
demonstrate how this might be done we first discuss the case M = 4. We then comment 
briefly on M > 4. 

We now define a 4-species model and argue that it phase separates. Consider a ring 
where each site is occupied by either an A, B, C or D particle. The model evolves according 
to the following procedure: at each time step two nearest neighbors are chosen randomly 
and exchanged according to the rates, 



As before the model conserves the number of particles of each species. Note that several 
other generalizations of the model to four species are possible. However, for simplicity, we 
discuss only the model defined by (|63"D with q < 1. 

We now argue that the system phase separates into a configuration of the form ABCD 
(where each letter now indicates a domain) as long as the densities of particles of each 
species are non-zero. Note that in the model AB, BC, CD, DA, AC and DB boundaries 
are stable while reverse boundaries BA, CB, DC, AD, CA and BD are unstable. As in 
the case of the three species model the system, starting from a random initial condition, 
evolves on a short time scale (i.e., which is not determined by the size of the system) into 
a metastable configuration where only stable domains are present. This configuration then 
slowly coarsens by slow diffusion of particles through neighboring domains. The system 
will finally reach the most stable state where the number of domains is minimal. One can 
easily check that this configuration is given by ABCD. Note that the system may exhibit 
other metastable states. For example, a state composed of ACDABCD is also stable under 
the choice of rates (|63|). However, since this state is composed of more domains than the 
4-domain state, some of the domains are necessarily smaller. According to the argument 
presented in Section II the relaxation time of this sequence (proportional to q~ m where m 
is the typical domain size) is much shorter than the relaxation time of the 4-domain state. 
Therefore the 4-domain state is more stable so that the system will finally evolve into it. 



VI. GENERALIZATION TO M SPECIES 




(63) 
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In considering M > 4 models one finds that for some choices of transition rates several 
states with the minimum number of domains may become metastable. For example, for 
M = 5 it is possible to choose transition rates for which both ABCDE and ACEBD are 
locally stable. The relative stability (and thus the resulting phase separated state) may be 
found by determining the relaxation time of these states using simple considerations such 
as those presented in Section II. 

As is the case of M = 3 detailed balance is found to be satisfied for certain densities and 
transition rates for M > 3. The condition for this is derived in Appendix B. In this case 
the relative stability of metastable states could be determined by comparing free energies. 



VII. CONCLUSION 

In this paper a model of three species of particles diffusing on a ring previously introduced 
in [T5|] has been studied. The model is governed by local dynamics in which all moves 
compatible with the conservation of the three densities are allowed. We argue that phase 
separation should occur as long as all densities are non-zero. In the special case of equal 
densities we find that the steady state generated by the local stochastic dynamics is exactly 
given by a long-range asymmetric Hamiltonian. Phase separation for this case is explicitly 
demonstrated. The model provides an explicit example for the mechanism leading to phase 
separation or breaking of ergodicity in systems with local stochastic dynamics. Although we 
did not succeed in solving the steady state in the case of non-equal densities of particles, 
there is a strong evidence that phase separation still occurs. In order to investigate further 
the case of equal densities we employed a generalized matrix ansatz to calculate correlation 
function for finite-size systems. The novel structure of the ansatz may give some clue as to 
handle other Id models which have so far resisted solution. 

The dynamics of phase separation reduces to a coarsening problem where the typical 
domain size grows logarithmically in time. This results from the elimination of the domains 
at a rate exponentially small in their size. The slow dynamics poses a problem of how to 
access the scaling regime numerically. With direct numerical simulations only small systems 
can be studied (see Fig. 3). However, by employing a toy model in which domains rather 
than individual sites are updated one can simulate much larger systems and probe the scaling 
regime (see Fig. 5). Such ideas of updating domains have been used before in the study of 



coarsening [pl]| . With the aid of the toy model it should be possible to study other aspects 
of the scaling regime associated with the slow dynamics and escape from metastable states. 

Generalizations are possible to models with M > 3 species. We have discussed some 
possibilities and have shown that phase separation may take place, although the structure 
of the set of metastable states is more complicated. As was the case for M = 3, conditions 
for detailed balance with respect to a long-range Hamiltonian may be determined. 

The problem of phase separation and coarsening is of interest also in the broader context 
of phase transitions in one dimensional systems. Here the existence of conserved quantities 
results in certain local transition rates being zero. It would be interesting to generalize this 
study to models in which no conserved quantity exists and all local rates are non- vanishing. 
Also, another open problem is to calculate the steady state of the present model in the case 
of non-equal densities. 
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APPENDIX A: MEAN-FIELD SOLUTION OF TOY MODEL 

To solve the mean-field model we follow the method used in PE| , PT|j . Let ni(t) be the 
number of domains of size I, irrespective of the type of particle it consists of. Let Z m in(£) 
be the length of the smallest domain and M(t) = 3K(t) be the total number of domains at 
time t. We will assume that ni(t) has the following scaling form in the large-time limit. 

M „, I , , k % 

ni = r-f(-r-) (Al) 

'min 'min 

After the elimination of the smallest domain as given by ([42"D , n h / min and M change according 
to 

M' = M- 3n Zmin (A2) 
n\ = n,(l - 5^p) + n lmin ^e(l - 2l min ) + 2n lmin '|f (A3) 

j 'min 

^min ^min 1 • (A4) 

Using the scaling form ( pTT| ) we have, 

< = 7^/(7-^) w J—[f( x ) ~ (3/(1) + l)y^ " j^d x f{x)\ (A5) 

'min t -L 'min T J- 'min 'min 'min 

where x = /// mm and we have expanded in l// m i n . Now substituting these in (|A3|) , it is 
straightforward to show that 



f(x) + xd x f(x) - 2/(x)/(l) + f(l)f(x - \)9{x - 2) + 26(x - 2)/(l) jf dy f(y)f(x - y) = . 

(A6) 

Using the Laplace transform, 

dx exp[— px] f(x) (A7) 

one can show that (f)(p) satisfies the differential equation, 

pd P (l>{p) = /(l)[0(p) - l][20(p) + e- p ] (A8) 



Since <p{p) = 1 — (x)p + . . ., by expanding (]A~8|) to order p one obtains /(l) = 1/3. The 
solution of (|A8|) with boundary conditions 0(0) = 1 and <p(p) « e~ p /3p for p >> 1 is 
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j'p doc x ii ^ 

W = 3p 2/3 e -/(p)/3 + dx piTi g _ x e _ /(a)/3 (A9) 

where I(x) = f£° dt exp[— xt]/t = — log(x) — 7 — J2™=i(~ x) n /{n\ n) and 7 is the Euler 
constant. Inverse Laplace transform of (|A9[) gives the domain size distribution. From ( |A9|) 
and using the expansion (p(p) = 1 — + . . ., where the average is with respect to f(x), 
we get 

3e 7/3 

^ = / °° dx x-V3 e- e-'W/ 3 ' (A10) 



APPENDIX B: DETAILED BALANCE CONDITION FOR AN M-SPECIES 

MODEL 

We now define the most general M species model, where M > 3. Let Xi denote a 
variable at site i of a ring of size iV, which takes values Xi = 1,2, ... , M. Xi = m means 
that site i is occupied by a particle of type m. The system evolves by a random sequential, 
nearest-neighbor exchange dynamics, with the following rates: 

q(m,n) 

mn ^2. nm , (Bl) 

q(n,m) 

and q(Xi,Xi) = 1. The model conserves N m , the number of particles of type m, for all m. 

We now present a condition for the model to satisfy detailed balance with respect to the 
steady- state weight given by 

JV-l N 

W({X l })= n II <l( X r- X i) : (B2) 

i=l j=i+l 

where the set {X^} describes the microscopic configuration. 

Consider a particle exchange between sites k and k + 1, where = m,Xfc +1 = n and 
k N (i.e. in the bulk, note that site 1 is chosen arbitrarily). Expanding the product in 
(|B2"D, it is easy to verify that 

W(Xi,...,m,n,...,X N ) = q(n,m) ^ 
W(Xx,...,n,m,...,X N ) q(m,n) ' 

Since this hold for any m, n, and is irrespective of the number of particles of each species, 
the steady-state weight ( |B2j ) satisfies detailed balance for all nearest-neighbor exchanges in 
the bulk. If the weights (|B2|) are translationally invariant then detailed balance will also 
hold for exchanges between sites 1 and N. 

Thus, to complete the proof of detailed balance it is sufficient to demand that (|B2| ) is 
translationally invariant. To do this we relabel sites i — > i + 1. The weight then becomes 

JV-l N 

w({x t })= n n (B4) 

i=l j=i+l 
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where X is identical to X^. Rewriting this equation by relabeling the indices we obtain, 



wax,}) 



N-l N 

n n viz*,* 

i=l j=i+l 



N-l 

n 



Comparing ([B5|) with ( p2|) and noting for example that, 



N M 

Y[q(X v X N ) = l[[q(l,X N )} 

j=l 1=1 

one can see that ( |B2| ) is translational invariant if 



A'; 



1=1 L 



q(l,m) 



(B5) 



(B6) 



(B7) 



for every m = 1, . . . , M. Thus, detailed balance holds if (|B7|) is satisfied. We note that for 
given densities {N m } the manifold of solutions for the rates is of M(M — 3)/2 dimensions. 
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